Graphene supports the propagation of subwavelength optical solitons 
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We study theoretically nonlinear propagation of light in a graphene monolayer. We show that 
the large intrinsic nonlinearity of graphene at optical frequencies enables the formation of quasi 
one-dimensional self-guided beams (spatial solitons) featuring subwavelength widths at moderate 
electric-field peak intensities. We also demonstrate a novel class of nonlinear self-confined modes 
resulting from the hybridization of surface plasmon polaritons with graphene optical solitons. 

PACS numbers: 42.65.Tg, 78.67.Wj, 73.20.Mf 



The experimental discovery and isolation of graphene 
monolayers from bulk graphite pQ has attracted great 
interest during the last years. The study of graphene 
properties has become a hot topic of research within the 
physics and nanoscience communities [2 as it promises, 
among others, a variety of optical and opto-electronical 
applications [3] |4] . Very large values of the nonlinear op- 
tical susceptibilities corresponding to multiple harmonic 
generation were theoretically predicted [5j [6] and have 
been experimentally verified very recently in the case 
of third-order nonlinear effects [7 . Still, it is an open 
question whether this high nonlinear coefficient, which 
occurs in a two-dimensional (2D) system, could induce 
strong nonlinear effects in electromagnetic (EM) modes 
that extend on the three spatial dimensions. 

One of the nonlinear effects with greater potential for 
controlling light propagation at the micro- and nano- 
scales is the formation of temporal and spatial EM 
solitons [8-12 . In this Letter we demonstrate that 
2D graphene monolayers support spatial non-diffracted 
beams (i.e., solitons) of subwavelength width in the opti- 
cal regime. We illustrate this capability by analyzing two 
arrangements leading to solitons with different polariza- 
tions: a graphene monolayer embedded into a conven- 
tional dielectric waveguide and a graphene sheet placed 
on top of a metal-dielectric structure. We analyze in 
detail the formation of spatial solitons and the relation 
between soliton width and input power, showing that 
the subwavelength scale can be reached by using feasi- 
ble values for the beam peak intensity. We also develop 
a quasi-analytical model that is able to capture the basic 
ingredients of the numerical results. 

The first structure in our analysis consists of a single 
graphene monolayer placed inside a planar linear dielec- 
tric waveguide, see Fig. [I] (a). This dielectric waveg- 
uide provides vertical confinement in the x-direction 
for the propagating EM mode. Graphene must be 
physically considered as a 2D material with nonlinear 
conductivity. But mathematically we can approximate 
graphene by a very thin layer of a finite thickness in- 




FIG. 1. (color online) Geometry and TE soliton forma- 
tion. An optical beam propagates inside the waveguide with 
a graphene monolayer located in the center in the low power 
(a) and high power (b) regimes. Panels (a) and (b) show 
slices of the beam intensity evaluated at the graphene layer, 
the yellow lines represent magnetic vector field whereas white 
lines depict the electric vector field. 



troducing an effective dielectric constant. Then we can 
treat graphene using the Maxwell equations for bulk me- 
dia. We have checked that both approaches give vir- 
tually the same numerical results. Since 2D and 3D 
treatments are equivalent, we take directly the nonlin- 
ear susceptibility from the experiment, and approximate 
graphene by a thin d gr = 0.3 nm- thick layer. The non- 
linear polarization density in graphene is Pnl(f, i) = 
fo^NiXr, £)Enl(Fj i), where Enl(f,£) is the electric field 
and 6NL(r,t) = Xg^r [Enl(f> t)] 2 is the nonlinear contri- 
bution to the equivalent permittivity of the graphene 
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layer. The parameter eo is the vacuum permittiv- 
ity. From the polarization density the nonlinear cur- 
rent is obtained as jNL(r,t) = SPnl(f, t)/dt. Through- 
out this paper, we consider the operating wavelength, 
Ao = 850 nm, for which a Kerr-type third-order ef- 
fective nonlinear susceptibility Xgr — 1.5 x 10 -7 esu, 
(xfj ~ 2.095 x 10~ 15 m 2 /V 2 in SI units) has been mea- 
sured [7]. The nonlinear eigenmode problem is formu- 
lated in terms of the 3D vector Maxwell equations and 
solved self-consistently in the continuous wave regime us- 
ing the finite element method [T3] . 

For the case depicted in Fig. 1(a), the initial solution 
for the iterative method has a form of a TE-polarized 
beam propagating in the z-direction with a gaussian 
shape along the ^/-direction and a waveguide profile in 
the x-direction. On each step of the iterative process, the 
EM fields are calculated by solving the propagation prob- 
lem with the eigenmode solution introduced as a source. 
In these calculations, we neglect third-order nonlinear ef- 
fects in the high-index dielectric material surrounding the 
graphene layer. This assumption is justified by the fact 
that the magnitude of the third-order nonlinear optical 
susceptibility in conventional high-index dielectric media 
is several orders of magnitude smaller than the one char- 
acterizing graphene. 

Figure 1 illustrates the formation of a non-diffracted 
beam at an operating wavelength Ao = 850 nm and 
for a dielectric waveguide of thickness 300 nm, char- 
acterized by a linear dielectric permittivity, = 2.25. 
When the beam intensity (defined as the modulus of the 
Poynting vector) at maximum is low (I < 10 13 W/m 2 ), 
the system operates in the linear regime and the beam 
diffracts while traveling in this structure, see Fig. 1(a). 
However, our numerical calculations show that, for high 
enough intensity (/ > 10 13 W/m 2 ), the nonlinearity of 
graphene can compensate diffraction leading to the for- 
mation of EM solitons. This is illustrated in Fig. 1(b), 
computed for / = 10 14 W/m 2 , showing a non-diffracted 
beam with a lateral size of the order of Ao- The soliton 
field is laterally confined due to the self-induced change 
of the effective refractive index, similarly to what hap- 
pens to a beam traveling within a bulk nonlinear waveg- 
uide [10] . In contrast to a conventional nonlinear waveg- 
uide, in which the nonlinear index change occurs in the 
whole volume, here in our system the 3D beam is lat- 
erally self-guided thanks to the nonlinearity that is only 
present in the 2D graphene sheet. We stress that these 
spatial solitons, sustained by a single graphene sheet, 
are very different to those supported by a metamaterial 
composed of graphene-dielectric superlattices in the ter- 
ahertz regime, which propagate perpendicularly to the 
graphene layers [14] . We also emphasize that the class of 
bright self-guided solitonic modes observed in Fig. 1(b) 
could not be supported by a thin metal film. In general, 
for the frequency range considered in this work, metal 
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FIG. 2. (color online) Soliton profile analysis. Panel (a) 
shows the transversal Enfield distribution, \E\, of the soliton 
mode for the case in which the peak intensity is / = 1.8 x 
10 14 W/m 2 . The horizontal cross-sections hi and /12 of the 
normalized Enfield (panel b) display a conventional soliton 
profile proportional to sech(y). For the vertical cross-sections 
(panel c) , the Enfield at V2 has the standard profile of a linear 
waveguide mode. For the cross-section evaluated at v\ the 
shape of the beam corresponds to that of a nonlinear system. 
The filled area on the panel (c) shows the location of the 
dielectric waveguide. 



films of nanometric thickness feature complex values of 
the nonlinear third-order susceptibility, x^ 3 \ such that 
Re^( 3 ) < and Imx^ is positive and large [15] . i.e., they 
display self-defocusing nonlinearities with large nonlinear 
absorption losses. 

Optical solitons in graphene should be observable with 
current samples and moderate beam intensities. Al- 
though the considered peak intensity in Fig. 1(b) is much 
higher than the reported damage threshold of graphene 
for continuous wave excitation, icwth ~ 10 10 W/m 2 [16], 
it is still well below the damage threshold of graphene for 
200 fs pulses, /pth ~ 10 16 W/m 2 [17 . Our continuous 
wave description of the soliton propagation under pulsed 
excitation is fully justified as the optical cycle associated 
with Ao = 850 nm is two orders of magnitude shorter 
than a 200 fs pulsed beam. 

Additional insight on the physical process can be 
gained from the study of the beam profile. Figure 2(a) 
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renders the normalized transversal electric field profile 
|E(#, of the calculated 3D eigenmode for / = 1.8 xlO 14 
W/m 2 (high power regime). The field cross-section along 
the ^/-direction [see Fig. 2(b)] can be accurately fitted 
by the function sech(y /w(x)), where w(x) is a measure 
of the lateral beam size, which slightly depends on x. 
This functional form for graphene EM solitons will be 
discussed later on. On the other hand, the confinement 
of the Enfield along the x-direction is governed by the to- 
tal internal reflection at the boundaries of the dielectric 
waveguide. The normalized E- field cross-section along 
the x-direction changes as the y coordinate is varied (see 
Fig.^c)). This change is more significant than in planar 
nonlinear waveguides and represents a distinct manifes- 
tation of the unusually large nonlinear optical current 
supported by the 2D graphene sheet that, in the present 
case, substantially exceeds the linear one. 

We turn now to analyze the dependence of soliton 
width (characterized by the full width at half maximum 
(FWHM), a, of the soliton E'-field) on the external in- 
tensity illuminating the system. In order to do this, 
we have computed the nonlinear eigenmodes for several 
peak E-^eld amplitudes in the graphene layer. The re- 
sults, in terms of the corresponding intensity distribu- 
tions (which in each case have been normalized to the 
maximum beam intensity), are summarized in the in- 
set of Fig. 3. As the maximum value of the E-field in 
the graphene layer (\E\ max ) is increased from 0.8 x 10 8 
V/m (bottom panel) to 4.2 x 10 8 V/m (top panel), a 
decreases from a=2 /im (more than 2 times the wave- 
length of the external illumination) to a=0.253 /im (well 
inside the subwavelength regime). The results displayed 
in Fig. 3 represent a novel instance, in a strict 2D system, 
on how the balance between nonlinearity and diffraction 
can yield self-guided propagating beams with subwave- 
length lateral confinement. In this context, it is impor- 
tant to point out that when graphene losses are incor- 
porated into the calculations (these losses stem from the 
linear part of the graphene conductivity), the propaga- 
tion length L of the soliton, defined as L = 1/2Iiii(/3/vl)> 
@nl being the complex propagation constant of the non- 
linear mode, is barely dependent on a. In fact, we have 
found numerically that L varies between 15 and 20 /im 
for all the soliton widths considered in this work. This 
independence of the propagation length on the field con- 
finement is very different to what is observed in other 
subwavelength-confined EM modes as, for example, sur- 
face plasmon polaritons. 

To account for the physical origin of the above de- 
scribed dependence of the soliton width on the peak elec- 
tric field amplitude, we have adapted to this problem the 
theoretical approaches used to describe soliton forma- 
tion in conventional 3D nonlinear optical materials [8]- 
[T0] . For this quasi- analytical treatment, we employ the 
3D modeling of the graphene layer, which, as mentioned 
before, gives virtually the same results as a description 
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FIG. 3. (color online) Dependence of the soliton width 
with the input power. The inset shows the normalized in- 
tensity distributions for decreasing values of the input power 
(from top to bottom) , resulting in spatial solitons of increased 
width. The operating wavelength in all cases is Ao = 850 nm. 
Main panel presents the peak electric field dependence with 
the soliton width. Circular dots represent the values obtained 
with the full nonlinear calculation whereas the solid line ren- 
ders the results from the quasi-analytical treatment (see main 
text). When the profile A(x) and the propagation constant 
f3 corresponding to the numerical calculation is included (tri- 
angular dots), the agreement between analytics and numerics 
is improved. 



based on a strictly 2D conductivity. Within this ap- 
proach the propagation of light inside the graphene layer 
is formulated in terms of the non-homogeneous vector 
Helmholtz's equation, 



c 2 e 



2 Q2 



A(r,*)=jtf L (r,t) (1) 



where A(r, t) is the magnetic potential vector (i.e., 
A(r,t) = -<9E(r, choosing the gauge V • A = 0) 

and n s is the linear refractive index of graphene. To 
solve Eq. (1), we start by assuming that its solutions are 
of the form 

A(r, t) = \ [A(x) F(z, y) exp[z(/3* - ut)] + c.c] (2) 



where A(x) is, in principle, an arbitrary function that 
governs the confinement of the EM field along the in- 
direction (see definition of axes in Fig. 1). As deduced 
from Eq. (2), A(x) also defines the polarization of the 
considered modal profile. The EM field profile in the 
graphene plane is controlled by the complex function 
F(z, y), whereas the corresponding propagation constant 
along the z-direction is given by /3. 

Now, we insert Eq. (2) into Eq. (1), we apply the 
slowly varying amplitude approximation, and we project 
the left and right-hand-side of the resulting equation 
over [A(x)]* T (where [ ]* T stands for the transpose 
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conjugate). Then, we define the auxiliary function 
f(z,y)=F(z,y)exp(-i</>z) (where <j> = (k 2 - 1 + 
J 2 //i)/2, fc s = n^ 2 /c 2 , h = fZodx\A(x)\ 2 and J 2 = 
J^oq dx[A(x)]* T d 2 A(x)/dx 2 ). Using these definitions, 
after some algebra, one finds that Eq. (1) can be rewrit- 
ten in terms of the function /(z, y) as 

^^S^^ ^I^C-, ^)|-^"C-, ^) = (3) 

where g = \u A X f) h/ he 2 , with J 3 = /^;/ 2 2 ^|A(x)| 4 . 
The crucial point to realize is that Eq. (3) corresponds to 
the standard form of the nonlinear Schrodinger equation, 
whose solutions have a canonical first-order soliton form 

f(y, z) = — J- sech(y/w) exp(iz/2f3w 2 ) (4) 
w\l g 

where w is the conventional definition of the soliton 
width, which in terms of the soliton FWHM is given 
by w ~ a/2.64. Physically, Eq. (3) and its correspond- 
ing solution given in Eq. (4) can be interpreted as those 
governing the propagation of light in a special class of 
index-guided waveguide in which the refractive index 
contrast between the core and the cladding is induced 
by the intensity of the propagating beam itself. Impor- 
tantly, Eq. (4) confirms the existence of soliton solutions 
in graphene, as observed in the numerical experiments 
reported in Figs. (l)-(3). Notice that the strength of the 
effective nonlinearity is characterized by the parameter g, 
which is proportional to both \gr (related to the intrin- 
sic nonlinearity of graphene in a free-standing configura- 
tion) and I3/I1 , which provides a measure of the fraction 
of EM energy that flows inside the graphene sheet. 

Inspired by the theoretical approaches used tradition- 
ally in nonlinear optics [SHIP], we assume that both the 
vector function A(x) and the propagation constant j3 of 
the modal profile correspond to those obtained numeri- 
cally for the linear counterpart of the structure sketched 
in Fig. 1(a). The results computed within this approxi- 
mation are displayed in Fig. 3 (see solid line), showing a 
qualitative agreement between the analytical results and 
the full numerical calculations. We emphasize that no 
fitting parameters are used in this comparison. The dis- 
crepancy between analytics and full numerics becomes 
larger as the value of a decreases. This fact can be 
ascribed to the difference between the profile A(x) ob- 
tained for the linear case and that computed numerically 
for the full nonlinear problem, which increases as a de- 
creases. This point is confirmed by the additional results 
displayed in Fig. 3 (triangular points), which show how 
the agreement between analytics and numerics improves 
when we introduce in Eq. (2) both the self-consistent pro- 
file, A(x) (obtained at y = 0), and the propagation con- 
stant j3 corresponding to our nonlinear simulations. The 
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FIG. 4. (color online) TM soliton formation. The geom- 
etry consists of the graphene monolayer and gold half-space 
separated by a silicon dioxide layer of thickness of 100 nm. 
The operating wavelength is Ao = 850 nm. The density plot 
in panel (a) is the transversal E-field distribution, \E\, asso- 
ciated with the excitation of a TM SPP-soliton. Panel (b) 
shows the dependence of the soliton FWHM, a, with the in- 
put power, measured as the peak Afield amplitude evaluated 
at the graphene monolayer. The dots represent the numerical 
results whereas the solid line is a fitting to a 1/a function. 
Inset of panel (b) renders an horizontal cross-section of the 
Enfield amplitude along the ^/-direction at the graphene layer. 



remaining difference can be traced back to the separabil- 
ity in the x and (y, z) coordinates implied by Eq.(2) that 
cannot fully account for the complexity of the graphene 
EM solitons. 

Finally, we show that TM-polarized optical solitons 
can also propagate along a graphene monolayer. A 
graphene structure that is able to support these TM op- 
tical solitons is rendered on Fig. [4] Here the vertical 
confinement is provided by a surface plasmon polariton 
(SPP) mode that is propagating on the interface between 
gold and a dielectric film. The graphene monolayer, 
which is characterized by a large nonlinear third-order 
susceptibility, must be separated from the metal surface 
by a dielectric spacer. We have chosen a 100 nm silicon 
dioxide layer, just for proof-of-principles purposes. Our 
calculations show that this system supports the propaga- 
tion of a very peculiar class of TM soliton, which results 
from the hybridization between the SPP supported by 
the metal-dielectric interface and the soliton propagat- 
ing in the graphene sheet. The computed transversal 
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E'-field distribution, |E(x,t/)|, of this hybrid SPP-soliton 
solution is plotted in Fig. [4ja) and displays exactly the 
conventional solitonic profile along the ^-direction, see 
the inset of Fig. 4(b). The dependence of the soliton 
width with the peak E'-field amplitude rendered in Fig. 
4(b) is very similar to that found for TE optical solitons, 
predicting the existence of subwavelength optical solitons 
also for this polarization. 

In conclusion, we have demonstrated that graphene 
monolayers can support both TE and TM spatial op- 
tical solitons due to the extremely large magnitude of 
its nonlinear third-order susceptibility. Moreover, we 
have shown that for feasible values of the input power 
these quasi-one dimensional optical solitons can have a 
subwavelength lateral width. We have also developed a 
quasi-analytical model that has a semi-quantitative value 
and that is able to predict the field intensities needed for 
soliton formation. The existence of subwavelength opti- 
cal solitons adds a new capability to the already broad 
range of optical phenomena associated with graphene 
structures. 
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